Exercise week 14
Shawn Schneidereit
2021-07-18
Data wrangel
data(UCBadmit, package = "rethinking")
d <- UCBadmit %>%
mutate(gid = ifelse(applicant.gender == "male", "1", "2"))
rm(UCBadmit)
d## dept applicant.gender admit reject applications gid
## 1 A male 512 313 825 1
## 2 A female 89 19 108 2
## 3 B male 353 207 560 1
## 4 B female 17 8 25 2
## 5 C male 120 205 325 1
## 6 C female 202 391 593 2
## 7 D male 138 279 417 1
## 8 D female 131 244 375 2
## 9 E male 53 138 191 1
## 10 E female 94 299 393 2
## 11 F male 22 351 373 1
## 12 F female 24 317 341 2# models
b11.7 <-
brm(data = d,
family = binomial,
admit | trials(applications) ~ 0 + gid,
prior(normal(0, 1.5), class = b),
iter = 2000, warmup = 1000, cores = 4, chains = 4,
seed = 42)
b11.8 <-
brm(data = d,
family = binomial,
bf(admit | trials(applications) ~ a + d,
a ~ 0 + gid,
d ~ 0 + dept,
nl = TRUE),
prior = c(prior(normal(0, 1.5), nlpar = a),
prior(normal(0, 1.5), nlpar = d)),
iter = 4000, warmup = 1000, cores = 4, chains = 4,
seed = 42)
# Creating custom distribution for m12.1
beta_binomial2 <- custom_family(
"beta_binomial2", dpars = c("mu", "phi"),
links = c("logit", "log"), lb = c(NA, 2),
type = "int", vars = "vint1[n]"
)
stan_funs <- "
real beta_binomial2_lpmf(int y, real mu, real phi, int T) {
return beta_binomial_lpmf(y | T, mu * phi, (1 - mu) * phi);
}
int beta_binomial2_rng(real mu, real phi, int T) {
return beta_binomial_rng(T, mu * phi, (1 - mu) * phi);
}"
stanvars <- stanvar(scode = stan_funs, block = "functions")
b12.1 <-
brm(data = d,
family = beta_binomial2, # the custom likelihood
admit | vint(applications) ~ 0 + gid,
prior = c(prior(normal(0, 1.5), class = b),
prior(exponential(1), class = phi)),
iter = 2000, warmup = 1000, cores = 4, chains = 4,
stanvars = stanvars, # note our `stanvars`
seed = 42)
b11.7_dept <- brm(data = d,
family = binomial,
admit | trials(applications) ~ 0 + gid + (0 + gid | dept),
prior(normal(0, 1.5), class = b),
prior(exponential(1), class = sd, group = dept),
prior(lkj(2), class = cor, group = dept),
iter = 2000, warmup = 1000, cores = 4, chains = 4,
seed = 42) models <- list(b11.7, b11.8, b12.1, b11.7_dept)
pp_models <- list(b11.7, b11.8, b11.7_dept)
purrr::map(pp_models, pp_check)## Using 10 posterior samples for ppc type 'dens_overlay' by default.
## Using 10 posterior samples for ppc type 'dens_overlay' by default.
## Using 10 posterior samples for ppc type 'dens_overlay' by default.## [[1]]
##
## [[2]]
##
## [[3]]
purrr::map(models, plot)
## [[1]]
## [[1]][[1]]
##
##
## [[2]]
## [[2]][[1]]
##
## [[2]][[2]]
##
##
## [[3]]
## [[3]][[1]]
##
##
## [[4]]
## [[4]][[1]]
purrr::map(models, base::summary)## [[1]]
## Family: binomial
## Links: mu = logit
## Formula: admit | trials(applications) ~ 0 + gid
## Data: d (Number of observations: 12)
## Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
## total post-warmup samples = 4000
##
## Population-Level Effects:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## gid1 -0.22 0.04 -0.30 -0.15 1.00 3753 2680
## gid2 -0.83 0.05 -0.93 -0.73 1.00 2663 2449
##
## Samples were drawn using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
##
## [[2]]
## Family: binomial
## Links: mu = logit
## Formula: admit | trials(applications) ~ a + d
## a ~ 0 + gid
## d ~ 0 + dept
## Data: d (Number of observations: 12)
## Samples: 4 chains, each with iter = 4000; warmup = 1000; thin = 1;
## total post-warmup samples = 12000
##
## Population-Level Effects:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## a_gid1 -0.53 0.53 -1.61 0.54 1.00 982 1196
## a_gid2 -0.43 0.54 -1.51 0.65 1.00 986 1238
## d_deptA 1.11 0.54 0.02 2.20 1.00 990 1232
## d_deptB 1.06 0.54 -0.02 2.15 1.00 995 1226
## d_deptC -0.15 0.54 -1.23 0.94 1.00 987 1266
## d_deptD -0.19 0.54 -1.26 0.91 1.00 991 1210
## d_deptE -0.63 0.54 -1.71 0.46 1.00 1004 1267
## d_deptF -2.19 0.55 -3.29 -1.08 1.00 1023 1336
##
## Samples were drawn using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
##
## [[3]]
## Family: beta_binomial2
## Links: mu = logit; phi = identity
## Formula: admit | vint(applications) ~ 0 + gid
## Data: d (Number of observations: 12)
## Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
## total post-warmup samples = 4000
##
## Population-Level Effects:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## gid1 -0.44 0.42 -1.25 0.37 1.00 3009 2283
## gid2 -0.32 0.44 -1.18 0.55 1.00 3301 2277
##
## Family Specific Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## phi 3.01 0.78 2.04 4.94 1.00 1909 1153
##
## Samples were drawn using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
##
## [[4]]
## Family: binomial
## Links: mu = logit
## Formula: admit | trials(applications) ~ 0 + gid + (0 + gid | dept)
## autocor ~ tructure(list(), class = "formula", .Environment = <environment>)
## Data: d (Number of observations: 12)
## Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
## total post-warmup samples = 4000
##
## Group-Level Effects:
## ~dept (Number of levels: 6)
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(gid1) 1.28 0.41 0.72 2.30 1.00 1279 1652
## sd(gid2) 1.55 0.49 0.88 2.77 1.00 1332 1735
## cor(gid1,gid2) 0.85 0.18 0.33 0.99 1.00 1102 2196
##
## Population-Level Effects:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## gid1 -0.52 0.47 -1.46 0.49 1.00 1140 1761
## gid2 -0.32 0.57 -1.44 0.82 1.00 1404 2112
##
## Samples were drawn using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).purrr::map(models, mcmc_plot)## [[1]]
##
## [[2]]
##
## [[3]]
##
## [[4]]
expose_functions(b12.1, vectorize = TRUE)## Trying to compile a simple C file## Running /Library/Frameworks/R.framework/Resources/bin/R CMD SHLIB foo.c
## clang -mmacosx-version-min=10.13 -I"/Library/Frameworks/R.framework/Resources/include" -DNDEBUG -I"/Users/shawn/Library/R/4.0/library/Rcpp/include/" -I"/Users/shawn/Library/R/4.0/library/RcppEigen/include/" -I"/Users/shawn/Library/R/4.0/library/RcppEigen/include/unsupported" -I"/Users/shawn/Library/R/4.0/library/BH/include" -I"/Users/shawn/Library/R/4.0/library/StanHeaders/include/src/" -I"/Users/shawn/Library/R/4.0/library/StanHeaders/include/" -I"/Users/shawn/Library/R/4.0/library/RcppParallel/include/" -I"/Users/shawn/Library/R/4.0/library/rstan/include" -DEIGEN_NO_DEBUG -DBOOST_DISABLE_ASSERTS -DBOOST_PENDING_INTEGER_LOG2_HPP -DSTAN_THREADS -DBOOST_NO_AUTO_PTR -include '/Users/shawn/Library/R/4.0/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp' -D_REENTRANT -DRCPP_PARALLEL_USE_TBB=1 -I/usr/local/include -fPIC -Wall -g -O2 -c foo.c -o foo.o
## In file included from <built-in>:1:
## In file included from /Users/shawn/Library/R/4.0/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp:13:
## In file included from /Users/shawn/Library/R/4.0/library/RcppEigen/include/Eigen/Dense:1:
## In file included from /Users/shawn/Library/R/4.0/library/RcppEigen/include/Eigen/Core:88:
## /Users/shawn/Library/R/4.0/library/RcppEigen/include/Eigen/src/Core/util/Macros.h:628:1: error: unknown type name 'namespace'
## namespace Eigen {
## ^
## /Users/shawn/Library/R/4.0/library/RcppEigen/include/Eigen/src/Core/util/Macros.h:628:16: error: expected ';' after top level declarator
## namespace Eigen {
## ^
## ;
## In file included from <built-in>:1:
## In file included from /Users/shawn/Library/R/4.0/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp:13:
## In file included from /Users/shawn/Library/R/4.0/library/RcppEigen/include/Eigen/Dense:1:
## /Users/shawn/Library/R/4.0/library/RcppEigen/include/Eigen/Core:96:10: fatal error: 'complex' file not found
## #include <complex>
## ^~~~~~~~~
## 3 errors generated.
## make: *** [foo.o] Error 1# "Unfortunately, using a custom distribution means that we can’t easily let STAN compute the ELPD for us, so we have to do it manually in R"
log_lik_beta_binomial2 <- function(i, draws) {
mu <- draws$dpars$mu[, i]
phi <- draws$dpars$phi
trials <- draws$data$vint1[i]
y <- draws$data$Y[i]
beta_binomial2_lpmf(y, mu, phi, trials)
}
# loo compair of all models
loo(b11.7, b11.8, b12.1, b11.7_dept)## Warning: Found 6 observations with a pareto_k > 0.7 in model 'b11.7'. It is
## recommended to set 'moment_match = TRUE' in order to perform moment matching for
## problematic observations.## Warning: Found 5 observations with a pareto_k > 0.7 in model 'b11.8'. It is
## recommended to set 'moment_match = TRUE' in order to perform moment matching for
## problematic observations.## Warning: Found 11 observations with a pareto_k > 0.7 in model 'b11.7_dept'. It
## is recommended to set 'moment_match = TRUE' in order to perform moment matching
## for problematic observations.## Warning: Not all models have the same y variable. ('yhash' attributes do not
## match)## Output of model 'b11.7':
##
## Computed from 4000 by 12 log-likelihood matrix
##
## Estimate SE
## elpd_loo -480.1 157.4
## p_loo 100.4 34.0
## looic 960.2 314.9
## ------
## Monte Carlo SE of elpd_loo is NA.
##
## Pareto k diagnostic values:
## Count Pct. Min. n_eff
## (-Inf, 0.5] (good) 3 25.0% 613
## (0.5, 0.7] (ok) 3 25.0% 172
## (0.7, 1] (bad) 1 8.3% 35
## (1, Inf) (very bad) 5 41.7% 1
## See help('pareto-k-diagnostic') for details.
##
## Output of model 'b11.8':
##
## Computed from 12000 by 12 log-likelihood matrix
##
## Estimate SE
## elpd_loo -57.3 8.6
## p_loo 12.4 3.6
## looic 114.6 17.1
## ------
## Monte Carlo SE of elpd_loo is NA.
##
## Pareto k diagnostic values:
## Count Pct. Min. n_eff
## (-Inf, 0.5] (good) 3 25.0% 2808
## (0.5, 0.7] (ok) 4 33.3% 267
## (0.7, 1] (bad) 3 25.0% 146
## (1, Inf) (very bad) 2 16.7% 18
## See help('pareto-k-diagnostic') for details.
##
## Output of model 'b12.1':
##
## Computed from 4000 by 12 log-likelihood matrix
##
## Estimate SE
## elpd_loo -68.9 3.0
## p_loo 2.0 0.6
## looic 137.8 6.0
## ------
## Monte Carlo SE of elpd_loo is 0.0.
##
## Pareto k diagnostic values:
## Count Pct. Min. n_eff
## (-Inf, 0.5] (good) 11 91.7% 1112
## (0.5, 0.7] (ok) 1 8.3% 990
## (0.7, 1] (bad) 0 0.0% <NA>
## (1, Inf) (very bad) 0 0.0% <NA>
##
## All Pareto k estimates are ok (k < 0.7).
## See help('pareto-k-diagnostic') for details.
##
## Output of model 'b11.7_dept':
##
## Computed from 4000 by 12 log-likelihood matrix
##
## Estimate SE
## elpd_loo -49.7 2.2
## p_loo 11.1 1.2
## looic 99.5 4.4
## ------
## Monte Carlo SE of elpd_loo is NA.
##
## Pareto k diagnostic values:
## Count Pct. Min. n_eff
## (-Inf, 0.5] (good) 0 0.0% <NA>
## (0.5, 0.7] (ok) 1 8.3% 299
## (0.7, 1] (bad) 6 50.0% 58
## (1, Inf) (very bad) 5 41.7% 11
## See help('pareto-k-diagnostic') for details.
##
## Model comparisons:
## elpd_diff se_diff
## b11.7_dept 0.0 0.0
## b11.8 -7.6 7.8
## b12.1 -19.2 1.4
## b11.7 -430.4 156.9Using the basic loo output we see that m11.7_dept has the lowest lepd. Interesting here is that its version b11.7 was the worst out of the 4 models, and by very large margin.
# previous method of waic and loo compare no longer seem to work due to the custom distribution of m12.1. The commented-out code is included below
waic_b11.7 <- add_criterion(b11.7, criterion = "waic")
waic_b11.8 <- add_criterion(b11.8, criterion = "waic")
waic_b12.1 <- add_criterion(b12.1, criterion = "waic")
waic_b11.7_dept <- add_criterion(b11.7_dept, criterion = "waic")
comparison_waic <- loo_compare(waic_b11.7, waic_b11.8, waic_b12.1, waic_b11.7_dept, criterion = "waic")
print(comparison_waic, simplify = FALSE)
loo_b11.7 <- add_criterion(b11.7, criterion = "loo")
loo_b11.8 <- add_criterion(b11.8, criterion = "loo")
expose_functions(b12.1, vectorize = TRUE)
# "Unfortunately, using a custom distribution means that we can’t easily let STAN compute the ELPD for us, so we have to do it manually in R"
log_lik_beta_binomial2 <- function(i, draws) {
mu <- draws$dpars$mu[, i]
phi <- draws$dpars$phi
trials <- draws$data$vint1[i]
y <- draws$data$Y[i]
beta_binomial2_lpmf(y, mu, phi, trials)
}
loo_b12.1 <- loo(b12.1)
loo_b12.1 <- add_criterion(b12.1, criterion = "loo")
loo_b11.7_dept <- add_criterion(b11.7_dept, criterion = "loo")
comparison_loo <- loo_compare(loo_b11.7, loo_b12.1, b12.1, loo_b11.7_dept, criterion = "loo")
print(comparison_loo, simplify = FALSE)